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The solid-solid coexistence of a polydisperse hard sphere system is studied by using the Monte 

Carlo simulation. The results show that for large enough polydispersity the solid-solid coexistence 
state is more stable than the single-phase solid. The two coexisting solids have different composi- 
tion distributions but the same crystal structure. Moreover, there is evidence that the solid-solid 
transition terminates in a critical point as in the case of the fluid-fluid transition. 
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The monodisperse hard sphere system is one of the best understood system in its equilib- 
rium phase behaviors, it undergoes an entropy-driven first-order transition from a disordered 
fluid to an ordered solid as the volume fraction increases |l|. The system often serves as 
an excellent starting point to study more complicated systems and a good model for the 
description of a class of colloidal dispersions. However, in a real colloidal system the par- 
ticles inevitably exhibit considerable size polydisperseity, which influences significantly the 
thermodynamic and dynamic behaviors of the system 2|, y|, |4|, |5|, |g, |7|, la, |9|, |l(j]. Therefore, 
a more realistic model describing hard sphere colloids is the size-polydisperse hard sphere 
system. The equilibrium phase behaviors of the polydisperse hard sphere system have not 
yet been fully understood. Besides the usual fluid-solid transition, a general consensus is 
that there exists a terminal polydispersity, above which the single-phase crystal becomes 
thermodynamically unstable. It is not clear what structure is the thermodynamically most 
stable when the polydispersity of the crystal exceeds the terminal polydispersity. Some 
theoretical studies show that beyond the terminal polydispersity the crystal will fractionate 



into two or more coexisting solid phases with the same crystal structure [J, la]. Others in- 
dicate that a crystahto-glass transition wi.l occur Q. Experimental.* the crystallisation 

does not occur at large enough polydispersity [2J. From the experiment we draw no defi- 
nite conclusion about the equilibrium phase behavior because the nonequilibrium effect will 
dominate the system for such a high polydispersity. So far, the solid-solid coexistence of 
polydisperse hard sphere system is only a theoretical prediction. In order to confirm the 
prediction, carefully designed experimental studies and comprehensive computer simulations 
are needed. To the best of our knowledge, the only simulation evidence up to now comes 
from the work by Fernandez et al. 12J. In that work they investigated the phase equilibria of 
the polydisperse soft-sphere system, and found that the crystal is highly inhomogeneous at 
large polydispersity. However, this does not provide us any details about the inhomogeneous 
structure. 

In this letter, we use the Monte Carlo method to investigate the solid-solid transition 
of a polydisperse hard sphere crystal. To simulate a polydisperse crystal we employ the 
semigrand ensemble, which is the best frame to study the polydisperse crystal [3|, lla, Il4l ]. 
In the ensemble the composition distribution is a priori unknown, and the independent 
variables are the chemical potential differences A/z(<r) of particles of each kind to a reference 
kind, which are given in advance. As a result, in the ensemble the conditions for coexistence 



of two phases will be satisfied if only the two phases have the same pressure and referenced 
chemical potential. As is well known, the first-order transition of a hard sphere system 
can be identified by looking for a double-peak structure in the volume histogram. Then, by 
tuning the pressure one can easily detect the transition point where the two peaks have equal 
weight [l5(. In actual computations we locate the transition point using a more tractable 



criteria "equal peak height" [16j, which differs from equal weight in the finite-size effect, and 
gives the same result in the thermodynamic limit. In the following we describe in detail the 
simulation method and discuss the obtained results. 

The quantity Pi SO (V) = Y(V)exp(— (3PV) measures the probability density to find a 
system with the volume V and the pressure P, here T(V) is the semigrand canonical partition 
function with given differences of excess chemical potential, defined as 

T(V) = NI Jn M S/"Sa Zn(V) X ^ | ^ ^ eM) " ^^ f II **' C 1 ) 

here, <Tj and a r are the diameter of the ith particle and an arbitrarily chosen refer- 
enced component, respectively. Zn is the canonical configuration integral Zn(V) = 
Iv " ' Iv e_/3C/ IL=i ^ r «' an d ^ex(o'i) = Ai(crj) — /cTln( — j^-) is the excess chemical potential 
relative to the ideal gas. We now consider the coexistence of two solids having the same 
structure, in this case the two solids can be connected with a reversible path as in the 
case of fluid-fluid coexistence. The path can be constructed in such a way by introducing 
an extended semigrand canonical ensemble with the volume as the ensemble variable. By 
changing the volume the system transforms from one coexisting solid to the other, provided 
that the solid-solid transition does exist. When we perform a random walk in the volume 
space, the probability density of finding the system in volume V is given by 

P ext (V)cxT(V,N,T,Afi ex (a)). (2) 

The semigrand canonical partition function T(V) can be simulated by using the flat his- 
togram methods 22j, |23|, [24| to within an overall multiplier, which is not needed in determine 
the relative probability. The simulation involves three kind of moves: particle displacement, 
particle resizing and volume changing. 

Previous calculations [IZj show that face-centered-cubic(fcc) phase is still the most sta- 
ble for the hard sphere crystal with a low size-polydispersity. Therefore, both coexisting 



solids are considered to be the fee structure. At present, we still have no definite knowl- 
edge about the position of polydispersity leading to instability of the single-phase crystal. 
Recent studies on the elasticity of the fee polydisperse hard sphere crystal found that there 
exists a mechanical terminal polydispersity(MTP) above which the crystal is mechanically 
unstable [1GJ. The MTP is an upper limit of the thermodynamical terminal polydispersity, 
thus we expect that a single-phase crystal is in the thermodynamically metastable state at 
polydispersity slightly below the MTP. The chemical potential difference function of the 
metastable crystal can be calculated with the SNERP algorithm [l7|, [l8|. And the solved 
chemical potential difference function is then used to study the solid-solid phase transition. 
To determine the MTP from elastic constants, a large amount of simulation time is needed. 
Here, we adopt a simple but effective approach, by which the MTP can be estimated roughly. 
The approach is based on the observation that, once the polydispersity is higher than the 
MTP, the SNERP procedure will converge very slowly or even does not converge. As a 
result, the MTP can be determined roughly by monitoring the chemical potential difference 
calculation with SNERP algorithm. Once we know the chemical potential difference func- 
tion of a metastable crystal, we can use it to calculate the probability densities P ext (V) and 
P±,o{V). 

We consider four initial metastable crystals at the same value for the composition distribu- 
tion, but with different volume fractions. Correspondingly, we have four different chemical 
potential difference functions. The initial prescribed composition is a truncated Schultz 
function, and the initial configuration is an ideal fee crystal. All the simulations are per- 
formed with a system of 256 size-polydisperse hard spheres in a parallelepiped box, periodic 
boundary conditions are used in all three directions. Figure [T](a) and figure [2j(a) are two 
typical chemical potential difference functions of the initial metastable crystals. They are 
obtained by applying the SNEPR method to the initial crystals. The initial crystal in Fig. 
[T]has higher volume fraction than the one in Fig. [2j so its chemical potential difference is 
larger. Figure [H(b) and figure [21(b) show the distributions of the dimensionless volume of the 
system, which are calculated from the extended semigrand canonical ensemble simulation by 
using the solved chemical potential difference functions plotted in Fig. Ufa) and Fig. [51(a), 
respectively. The dimensionless volume is defined as V* = V/N(a) 3 , here a is the average of 
the diameter of particles of the initial metastable crystal. The distributions clearly exhibit 
a double-peak structure for a range of pressure, which is the sign of the existence of solid- 













(a) 


^~^ 




8 




/ \ 


/ 


4 




/ \ 


jf 


i 

< 






- 


•4 

-8 


^^ 



0.004 


*\ I 


- 


0.003 




- 


> 0.002 
o 


- 


- 


D." 

0.001 


■' \ 


- 


0.000 


V I 


1: 



1.5 1.6 1.7 1.8 1.9 2.0 2.1 



0.70 0.77 0.84 0.91 0.98 1.05 1.12 

V* 



FIG. 1: (a) the chemical potential difference of an initial metastable crystal as function of the 
diameter of particles. It is obtained from the SNEPR method. The inset shows the prescribed 
composition distribution of the metastable crystal, (b) the distribution of the dimensionless volume 
of the system at the coexisting pressure. It is calculated from the extended ensemble simulation 
by using the A/j, ex (a) plotted in (a). 
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FIG. 2: The same as Fig. [H the volume fraction of the initial metastable crystal is less than that 
in Fig. [TJ 



solid coexistence. At the coexisting pressure the two peaks have equal height, each peak 
represents a coexisting solid, and the volumes of two coexisting solids are determined from 
the positions of the peak maximum. The composition distribution of each coexisting solid 
can be obtained from an additional semigrand canonical simulation by the same chemical 
potential difference function. Figure [3] displays the particle size distributions of the coex- 
isting solids. The compositions of two coexisting solids are significantly different, i.e. the 
fractionation effect is sufficiently large. The coexisting solid with lower volume fraction has 



05 


(a) 










f 




T1=0 


634 


6= 


=0.0334 


0.04 




I 






": 




0.03 


- n=0-578, 8=0.0503 


0.02 


^^J 










- 


0.01 
0.00 










/ 

















0.020 


(b) 






, 




0.015 


n= 


0.568 


5=0.0748 


0, 0.010 










/ ^^%^ 


0.005 








"1 


- 


0.000 


1= 


=0.594 


8=0.0573 



1.5 1.6 1.7 1.8 1.9 2.0 2.1 



1.5 1.6 1.7 1.8 1.9 2.0 2.1 



FIG. 3: (a) the composition distribution of the coexisting solids given in Fig. [2(b). The solids 
with the lower(open square) and higher(filled square) volume fractions correspond to the left and 
right peaks in Fig. QJb), respectively. This is because the large volume accommodates easily the 
larger particles. Here, r\ is the volume fraction and 5 is the polydispersity. (b) the composition 
distribution of the coexisting solids given in Fig. [2(b). 



a larger polydispersity than the one with higher volume fraction, which is consistent with 
the results by Fasolo et al [5|. In the simulation we find that the separation between the 
two peaks decreases as the volume fraction of initial solid decreases, as shown in Fig. QJb) 
and Fig. [2(b). We expect that the two peaks will completely merge together under some 
special conditions. This gives the solid-solid critical point, sketched by the filled circle in 
Fig. [U The early studies also indicated that a monodisperse system of hard spheres with a 
short-ranged attractive interaction undergoes a solid-solid transition with the same crystal 
structure, and has a solid-solid critical point 25|, [26|]. We argue that the critical phenomena 
of the polydisperse crystal is not experimentally observable, since in a real colloidal crystal 
the particles are not permitted to exchange their positions and to change their sizes. 

Figure H] shows the coexisting solids plotted in the polydispersity and volume fraction 
plane. Even though we do not provide the cloud points(The estimate of cloud points need 
more simulation time and techniques 27|.). from the figure we can draw some important 
conclusions. For the systems under consideration we note that, except the lower density 
solid in Fig. [3(b) (the uppermost triangle in Fig. 0J), the coexisting phases with lower 
volume fraction have similar size distribution function and follow a linear relation. So the 
linear fit roughly resembles a cloud curve(dashed line). On the right side of the cloud 
curve the single-phase solid becomes thermodynamically unstable. The cloud curve has a 
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FIG. 4: The solid-solid coexistence plotted in rj — 5 plane. Two coexisting phases are denoted 
by the same symbols. Dashed line is a linear fit for the coexisting solids with the lower volume 
fraction (except the uppermost triangle); dotted line is a linear fit for the coexisting solids with the 
higher volume fraction; filled circle gives a guess for the critical point. 



negative slope, so for each volume fraction of interest there exists a maximum polydispersity, 
above which the solid-solid coexistence occurs. On the other hand, all the coexisting solids 
with higher volume fraction also have the similar sizes distribution. Thus we get a second 
approximate cloud curve by a linear fitting(dotted line). Comparing with the first case, 
the unstable region lies on the left side of the cloud curve, and for each volume fraction 
under consideration there exists a minimum polydispersity stabilizing a single-phase solid. 
Because the monodisperse crystal is stable in the volume fraction range of the simulation, 
the cloud curve will bend back at the lower polydispersity, where the fluid-solid transition 
occurs. However, with further increasing the polydispersity the cloud curve will bend back 
toward the high volume fraction due to the existence of the terminal polydispersity. The 
cloud curve looks like a "reverse -Z". 

In summary, we investigated the solid-solid transition of the size polydisperse hard sphere 
system by the Monte Carlo simulation. The results indicate that at sufficient high polydis- 
persity the single-phase solid become unstable and the solid-solid coexistence will occur. 
Two coexisting solid phases have different composition distributions. We also find that the 
solid-solid critical point can exist under some special conditions. Though we do not give 
the cloud points, we can obtain some qualitative understanding about the solid-solid phase 
diagram. It is possible to calculate the cloud points in the current frame [27j. The existence 



of the solid-solid coexistence state can not preclude the occurrence of the glass transition, 
but we prefer to think of the glass transition as a kinetic phenomenon. 
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